// Copyright (C) 2009,2010,2011,2012  Marco Restelli
//
// This file is part of:
//   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
//
// LDGH is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// LDGH is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with LDGH. If not, see <http://www.gnu.org/licenses/>.
//
// author: Marco Restelli                   <marco.restelli@gmail.com>

// Mesh size
//h = 8.0; // implies 224 elements
//h = 4.0; // implies 906 elements
//h = 2.0; // implies 3780 elements
//h = 1.0; // implies 15118 elements
h = 0.5; // implies 59150 elements
//h = 0.25; // implies 238368 elements

h_cil = h/8;

// Define the cylinder
xcn = 0.0;
ycn = 0.0;
rad = 0.5;
// the cylinder is defined using four arcs
cen = newp; Point(cen) = {xcn,ycn,0.0,h};
A = newp; Point(A) = {xcn,ycn-rad,0.0,h_cil};
B = newp; Point(B) = {xcn+rad,ycn,0.0,h_cil};
C = newp; Point(C) = {xcn,ycn+rad,0.0,h_cil};
D = newp; Point(D) = {xcn-rad,ycn,0.0,h_cil};

a1 = newreg; Circle(a1) = {A,cen,B};
a2 = newreg; Circle(a2) = {B,cen,C};
a3 = newreg; Circle(a3) = {C,cen,D};
a4 = newreg; Circle(a4) = {D,cen,A};

hol = newreg; Line Loop(hol) = {-a4,-a3,-a2,-a1}; 

// Now the outer box
dxm =  5.0;
dxp = 20.0;
dy  =  7.0;
cbl = newp; Point(cbl) = {-dxm,-dy,0.0,h};
cbr = newp; Point(cbr) = { dxp,-dy,0.0,h};
ctr = newp; Point(ctr) = { dxp, dy,0.0,h};
ctl = newp; Point(ctl) = {-dxm, dy,0.0,h};

bot = newc; Line(bot) = {cbl,cbr}; 
rig = newc; Line(rig) = {cbr,ctr}; 
top = newc; Line(top) = {ctr,ctl}; 
lef = newc; Line(lef) = {ctl,cbl}; 

box = newreg; Line Loop(box) = {bot,rig,top,lef}; 

// Finally define the domain
domain = news; Plane Surface(domain) = {box,hol};

// Display boundary labels
View "Boundary labels" {
  T2(10, 20, 0){ StrCat( 
    "Hole labels:    ",Sprintf("%.0f, %.0f, %.0f, %.0f",a1,a2,a3,a4)) };
  T2(10, 40, 0){ StrCat( 
    "Inflow label:   ",Sprintf("%.0f",lef)) };
  T2(10, 60, 0){ StrCat( 
    "Outflow labels: ",Sprintf("%.0f, %.0f, %.0f",bot,rig,top)) };
};

// Central refinement
Field[1] = MathEval;
Field[1].F = Sprintf( "%g/2*(0.25+1.75*(y/%g)^2)" , h , dy );
Background Field = 1;

